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The fourth virial coefficient of asymmetric nonadditive binary mixtures of hard disks is computed with a 
standard Monte Carlo method. Wide ranges of size ratio (0.05 < q < 0.95) and nonadditivity (—0.5 < 
A < 0.5) are covered. A comparison is made between the numerical results and those that follow from some 
theoretical developments. The possible use of these data in the derivation of new equations of state for these 
mixtures is illustrated by considering a rescaled virial expansion truncated to fourth order. The numerical 
results obtained using this equation of state are compared with Monte Carlo simulation data in the case of a 
size ratio q = 0.7 and two nonadditivitics A = ±0.2. 



I. INTRODUCTION 

The key role that hard-core model systems play in 
liquid state theory is undeniable. This is mostly due 
to the well-known fact that in some cases it is possible 
to derive exact and approximate analytical results for 
their thermodynamic and structural properties^ More- 
over, the structural properties of real dense fluids depend 
essentially on the short ranged repulsive intermolecular 
forces, which are adequately accounted for by hard-core 
models in which molecules have no interactions at separa- 
tions larger than a given distance and experience infinite 
repulsion if their separation is less than that distance. 
While pure one-component hard-core systems lead to a 
fluid-solid transition, mixtures may display more com- 
plex phase behavior. For these latter, one can either 
assume that they are additive, namely that the closest 
distance of approach of molecules of two different species 
is the arithmetic mean of the distances between like pairs, 
or nonadditive, in which the previous condition does not 
hold. Additive systems have received most of the atten- 
tion, but the inclusion of nonadditivity, which may be 
either positive or negative, attempts to incorporate some 
features of non-hard forces, such as attractions and soft 
repulsions, into the description. Amongst other things, 
nonadditivity serves to account for homo-coordination or 
hetero-coordination in the compositional order of a mix- 
ture and also for fluid- fluid demixing. This makes the 
nonadditive hard-core models of mixtures both attrac- 
tive and rather versatile and so it is not surprising that 
they have been the subject of recent attention in the 
literature. Some examples concerning nonadditive hard 
spheres (NAHS) may be found in Refs. 
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As far as mixtures of nonadditive hard disks (NAHD) 
are concerned, which are the subject matter of this pa- 
per, publications are less numerous than in the case of 
NAHS. However, interest in these model systems, which 
dates back at least to the late 1970s, has recently ex- 
perienced a revival. Applications include lipid monolay- 
ers spread on air- water interfaces^ liquid-liquid demixing 
in a physisorbed mixture of Argon, Krypton, or Xenon 
on graphite^ a model for ganglioside lipid and phos- 
pholipid interactions in connection with the binding of 
cholera-toxin to a lipid membrane,— the morphology of 
composite latex particles^ two-dimensional magnetic col- 
loid mixtures^ and the asphaltene flocculation inhibi- 
tion phenomenon^ 

A binary mixture of NAHD is characterized by the im- 
penetrable diameters of the two species an = &\ and 
022 = 02 and by a crossed diameter a 12 = 5(01 + 
02) (1 + A), where the dimensionless parameter A ac- 
counts for deviations of the inter-species interactions 
from additivityJ^ Like in the NAHS model, the binary 
mixture shows a tendency to form hetero-coordinatcd 
clusters for negative values of the nonadditivity pa- 
rameter (A < 0). On the other hand, for positive 
non-additivity (A > 0), the system tends to segre- 
gate into two fluid phases, one richer in particles of 
species 1 and the other richer in particles of species 2, 
respectively^ On the computational side, Dickinsons- 
reported molecular dynamics simulations of NAHD mix- 
tures in which he computed the compressibility factor 
and the radial distribution functions for a few size ratios 
and some nonadditivities. Tenne and Bergmanoai devel- 
oped a scaled-particle theory (SPT) for NAHD mixtures 
which was later corrected by Bearman and Mazo^— 
in their study of fluid-fluid phase equilibria for positive 
nonadditivity. The compressibility factors and part of 
the coexistence curve arising from the SPT were com- 
pared to molecular dynamics simulations of an equimo- 
lar symmetric mixture of NAHD by Ehrenberg et al£^ 
Singh and Sinha^ 2 . used thermodynamic perturbation 
theory to compute the Helmholtz free energy per parti- 
cle, the compressibility factor, and the radial distribution 
function of binary NAHD mixtures with both positive 
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and negative nonadditivity, while Mishra and Sinha£2 
derived the excess thermodynamic properties of binary 
NAHD mixtures including quantum corrections. Niclaba 
and coworkers^— combined the Gibbs ensemble Monte 
Carlo (GEMC) method and finite-size scaling to study 
demixing of a symmetric NAHD mixture. Hamad and his 
collaborator ;) 28 ' 29 developed equations of state for NAHD 
mixtures and performed molecular dynamics simulations 
for a variety of size ratios and values of the nonadditiv- 
ity. Saija and Giaquinta— reported Monte Carlo (MC) 
results for the thermodynamic and structural properties 
of a symmetric NAHD mixture for positive nonadditivity 
and studied phase separation for some positive values of 
the nonadditivity. Depletion interactions in NAHD mix- 
tures were considered by Castaheda-Priego et aZ.,— who 
also indicated that this model may mimic the qualita- 
tive features of effective potentials of hard and soft par- 
ticles. To cope with large nonadditivitics, Buhot^I used a 
cluster algorithm to study phase separation of symmetric 
binary NAHD mixtures, while Guaqueta^ 2 - used a combi- 
nation of MC techniques to determine the location of the 
critical consolute point of asymmetric NAHD mixtures 
for a wide range of size ratios and values of the pos- 
itive nonadditivity. More recently, Muhoz-Salazar and 
Odriozola^ used a semi-grand canonical ensemble Monte 
Carlo method to obtain the fluid-fluid coexistence curve 
for a symmetric mixture of NAHD and a single positive 
nonadditivity. 

In 2005 three of us^i introduced an approximate equa- 
tion of state for nonadditive hard-core systems in d di- 
mensions and, taking d = 2, compared the results ob- 
tained for the corresponding compressibility factor with 
simulation data. Later, a unified framework for some of 
the most important theories (including some generaliza- 
tions) of the equation of state of d-dimcnsional nonaddi- 
tive hard-core mixtures was presented^ The framework 
was used for d = 3 to compare the results of the different 
approaches with simulation data for the fourth virial co- 
efficients that had recently been derived^ and with sim- 
ulation data for the compressibility factor. It was also 
used to examine the issue of fluid-fluid demixing. 

More recently, another of us^I computed the fourth 
virial coefficient of symmetric NAHD mixtures over a 
wide range of nonadditivity. He also compared the fluid- 
fluid coexistence curve derived from two equations of 
state built using the new virial coefficients with some 
simulation results. 

One of the major aims of this paper is to present the 
results of computations of the fourth virial coefficient of 
asymmetric NAHD mixtures, i.e., mixtures such that the 
size ratio q = Gijo\ is different from unity. We will ex- 
plore a wide range of values of the nonadditivity parame- 
ter A and size ratio q. These results complement the ones 
already published for symmetric mixtures 3 ^ and will af- 
terwards be used to assess the merits and limitations of 
some theoretical approaches. 

The paper is organized as follows. In section [TT] we 
provide the known analytical results for the second and 



third virial coefficients of a NAHD mixture, as well as 
the graphical representation of the (partial) composition- 
independent fourth virial coefficients. The approximate 
theoretical expressions considered in this paper for the 
fourth virial coefficients are presented in section lHTl This 
is followed in section IIVI by the results of the MC evalu- 
ation of the fourth virial coefficients for a wide range of 
size ratios and values of the nonadditivity parameter. A 
comparison of the theoretical approximations with these 
data is also presented. In section fVl the equation of state 
resulting from a rescaled virial expansion truncated to 
fourth order, as well as the theoretical approximations 
mentioned above, are compared with new Monte Carlo 
simulation data in the case of two mixtures with negative 
and positive nonadditivitics, respectively. The paper is 
closed in section IVT1 with some concluding remarks. 

II. VIRIAL COEFFICIENTS 

The virial expansion can be written as 

PP = p + Bp 2 + Cp 3 + Dp 4 + ■ ■ • , (1) 

P is the pressure, /3 is the inverse temperature in units of 
the Boltzmann constant and p = p± + p 2 is the total num- 
ber density, pi being the partial number density of species 
i. In a mixture, at variance with the one-component case, 
the virial coefficients B,C,D,. . . do also depend on the 
relative concentration of the two species and on the hard- 
core diameters. The coefficients B and C are exact and 
well known (see, for instance, Refs. l28l and l34h . They are 
given by 

B = B n x 2 + 2Bi 2 xix 2 + B 2 2xl, (2) 

C = Cl\ix\ + iC\\ix\x2 + 3Ci22Xlxl + C 2 22xl, (3) 

where x\ = pi/p and X2 = Pil P = 1 — %i are the mole 
fractions of species 1 and 2, respectively. The other quan- 
tities read 

Bij = l4' ( 4 ) 
du = ^b 3 af, (5) 

7T 2 

C222 = T^3C2i (8) 
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where 6 3 = ^ - ^ ~ 3.12802 and the function F(x) is 
given by 

F(x) = ^G(x) + ^x 2 H(x) (9) 

with 

G(x) = J- (ax 2 cos" 1 i- - v 7 ^ 2 - 1^ , (10) 



contributes a factor to the integrand in the form of a 
Mayer step function. Space integration is carried out over 
all the vertices of the graph. Of course, the coefficient 
-D1222 is obtained from Eq. (|16p by exchanging the open 
and solid circles. 

For later use, let gij (p) be the values of the radial distri- 
bution functions at contact of the NAHD mixture. This 
quantity is related to the pressure via the virial equation 
of stated 



7T0 3 



2ttx 2 - 2 (2x 2 - l) cos 



for x > g and 



for < x < \. 



G(x) = 0, H{x) = —x 2 
03 



- 1 
(11 



(12) 



In turn, the fourth-order virial coefficient reads 

D = Dmixf + ADm 2 x\ x 2 + §Dn 22 x\x% 
+AD 122 2Xixl + D 2222 x 2 , 



(13) 



and its partial contributions have to be evaluated nu- 
merically. The terms Dim and D 2222 can be calculated 
through the expression of the fourth virial coefficient for 
a monodispcrsc fluid of particles with diameter ai or a 2 , 
respectively, i.e., 



PP = P + JP 2 X i X J a ij9ij(P)- 



(18) 



No general expression is known for gij(p), but it may 
formally be expanded in a power series in density as 



9ij 



k=l 



7T * 7T 2 ^ 

(p) = 1 + - 1 >y2-n. <■>.:,, + y^p 2 Y x k x e d ke . 



k,l=l 



(19) 

where the coefficients c^ij , dkt.ij , • • ■ are independent of 
the mole fractions but in general depend in a non-trivial 
way on the set of diameters {or/}- Only the coefficients 
linear i n p (i.e., Ck-ij) are known analytically (c/. Refs. 
|28| and I34 ) , namely 
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Dim = g^^o - !, 



7T 

64' 



D 2222 = —b±(j 2 , 



(14) 



(15) 



where 64 = 8(2 + 10/tt 2 - 9V3/2tt) ~ 4.25785. On the 
other hand, the coefficients Dm 2 and D1122 are clus- 
ter integrals which are represented by the following four- 
point color graphs: 



Di 
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0- 



(16) 



(17) 



The open and solid circles in each graph identify particles 
belonging to species 1 and 2, respectively. Each bond 



b 3 2n (<?i2\ b 3 2 (<t 12 \ 

C2;11= ^W)' Cl;12 = -2 aiH {-^)- (21) 

Other combinations of indices follow from the exchange 
of indices 1 and 2 in the above results. We recall that the 
functions G(x) and H(x) are given by Eqs. (fT0|) — (fT2 |) . In 
fact, insertion of Eq. (fH?)) into Eq. (TT8"|) yields 



C W = ^ {°k;ij^ij + Cj-tk^k + c i;jkVjk) ) ( 22 ) 

so that Eqs. ©--((8]) are recovered from Eqs. ([20| and 



III. APPROXIMATE THEORETICAL APPROACHES 

Before we evaluate numerically the partial fourth virial 
coefficients, let us recall the approximate results derived 
for them with different theoretical approaches. These 
were presented in a unified framework within the descrip- 
tion of general multi-component nonadditive hard-sphere 
mixtures in d dimensions. We will consider here the par- 
ticular case of a binary mixture in two dimensions, only 
quote the relevant results, and refer the interested reader 
to Ref. K59 for details. 
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A. MIX1 approximation 



D. The Santos-Lopez de Haro-Yuste proposal 



In the so-called MIX1 theory for NAHD mixtures, 
which wc will label with a superscript M, the fourth virial 
coefficients are given by 

2 
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7add 

a ik;j( 



1 + 3K 



(23) 



In Eq. (j2"3"j) , are the second-order coefficients defined 
in Eq. (|19p. particularized to the additive case (A = 0). 
Here we adopt the approximatio n 3 



jadd 

due.;; = Oi.O 



2^2 



l kt;ij 

Moreover, in Eq 



2 / Pi + U; <J k (T£ 



r^ = 2A(i-<y tf ), 



(24) 



(25) 



where 5{j is the Kronecker delta. 



B. Paricaud's modified MIX1 theory (mMIXl) 

In the generalization of Paricaud's approximation that 
was made in Ref. [HI, which will be identified with 
the superscript mM, and restricting the result to two- 
dimensional binary mixtures, the partial composition- 
independent fourth virial coefficients have the same form 
as in the MIXf approximation but one has to replace Y^ 1 
with >7 nM , where this latter is given by 



r- M ^A(2 + A)(I-^). 



(26) 



C. Hamad's proposal 

In the work of Hamad and his collaborators , 28 ' 29 de- 
noted here by the superscript H, the fourth virial coeffi- 
cients are given by 

t~~iH ^4 ( 2 2 2 



2 2 2 \ 

+ a jk c i:,jk c £:jk + & jt C i;jt c k;jl + &ki c i;klCj;kt) 



(27) 



In the proposal made in 2005 by three of us,— here- 
after denoted by the superscript SHY, the fourth virial 
coefficients arc expressed in terms of the partial second 
and third composition-independent virial coefficients and 
of 63 and 64. Written for d = 2 they read 



D 



SHY _ MPi 2) f2(j.. , 

16(6 3 - 2" 



ijk£ 



jki 



k^ijt 



7T 2 (6 4 - h) , 2 2 . 2 2 R ,2 2 R 



96(6 3 - 2) 
+a 2 ;a\Bu + a^crfBik + a\a\B l 



(28) 



IV. RESULTS 



In this section we report the results of our calculations. 
In order to evaluate the irreducible cluster integrals which 
enter the expression of the composition-independent co- 
efficients Dijki [see Eqs. (flU)) and (fTT)) ]. we used a stan- 
dard MC integration procedure. The algorithm produces 
a significant set of configurations which are compatible 
with the Mayer graph one wants to evaluate. We first fix 
particle 1 of species i at the origin and sequentially de- 
posit the remaining three particles at random but in such 
a way that particle a + 1 overlaps with particle a (where 
a = 0,1,2,3). This procedure generates an open chain 
of overlapping particles which is taken as a "trial config- 
uration" . A "successful configuration" is a closed-chain 
configuration (i.e., a configuration in which particle 1 fur- 
ther overlaps with particle 4) where, moreover, the resid- 
ual cross-linked "bonds" which are present in the Mayer 
graph that is being calculated are also retrieved. The ra- 
tio of the number of successful configurations (N s ) to the 
total number of trial configurations (N t ) yields asymp- 
totically the value of the cluster integral relative to that 
of the open-chain graph which, in turn, is trivially related 
to a product of the partial second-order virial coefficients 
Bij i 42 ' 43 The numerical accuracy of the MC results obvi- 
ously depends on the total number of trial configurations. 
The error on the cluster integral J is estimated as4i 



error = 



J(J-l) 



-1 1/2 



Ni 



(29) 



However, as a result of the accumulation of statistically 
independent errors, the global uncertainty affecting the 
partial virial coefficients is higher than the error esti- 
mated for each cluster integral that enters the expression 
of Dij k £. A typical MC run consisted of 4 x 10 9 inde- 
pendent moves. The error on each cluster integral, as 
estimated through Eq. (|2"9")l , turned out to be systemati- 
cally less than 0.05%, with a cumulative uncertainty on 
the partial virial coefficients lower than 0.5%. 

The numerical values of -Dn^/ci, -D1122/011 an d 
£1222/01 for A = ±0-05, ±0.1, ±0.2, ±0.3, ±0,4, 
and ±0.5 and q = 0.05, 0.10, . . . , 0.90, 0.95 are presented 
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FIG. 1. Plot of the composition-independent fourth virial 
coefficients D1112, -D1122 and -D1222 versus the size ratio 
q — 02/cri for a nonadditivity parameter A = —0.1. The 
dotted (green) lines correspond to the original MIX1 theory, 
Eq. (|23|) , the dash-dot (pink) lines correspond to the mMIXl 
theory, Eq. (g3J, with Y$ -> Y™ M , the dashed (blue) lines 
correspond to Hamad's proposal, Eq. (|27|l , and the solid (red) 
lines correspond to the SHY proposal, Eq. (|28|l . The symbols 
are our MC data. 



in tabular form in the supplementary material to this 

45 

paper. 

Now we proceed to assess the merits of the differ- 
ent theoretical formulae for the composition-independent 
partial fourth virial coefficients that we presented in sec- 
tion Mil For that purpose, although we have made an 
exhaustive analysis, in Figs.[T]-[S]we present only some il- 
lustrative cases in which we compare the performance of 
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cx/a, 

FIG. 2. Same as in Fig. [1] but for A = -0.3. 



the different approximations against the MC data. The 
graphs corresponding to the other values of A that ap- 
pear in the tables of the supplementary material to this 
paper— are available upon request. 

From these figures it is clear that, overall, the proposal 
by Hamad^ Eq. (|27|) . is very good for -D1112 and -D1222 
but rather bad for -D1122 if |A| > 0.1, irrespective of the 
value (positive or negative) of A. None of the theories 
shows a good performance in the case of -D1122 but at 
least the SHY proposal leads to reasonable quantitative 
agreement in the positive region of this coefficient, be- 
ing particularly superior to all other approximations for 
negative values of A. 
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FIG. 3. Same as in Fig. [TJ but for A = -0.5. 
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FIG. 4. Same as in Fig. [2 but for A = 0.1. 



fori 



V. EQUATION OF STATE 



Since the convergence of the virial expansion is un- 
known and truncating the series after the first four terms 
would not guarantee a satisfactory outcome, in this sec- 
tion we will use the knowledge of the first four virial 
coefficients to illustrate the performance of a well estab- 
lished approach to the equation of state of fluids that 
incorporates such knowledge. Hence we will consider the 
rescaled virial expansion (RYE) proposed by Baus and 
Colo t 46 i 47 to obtain an (approximate) equation of state 
for an asymmetric NAHD mixture. The RVE equation 
of state truncated to the fourth order has the following 



I3P _ 1 + ci?7 + c 2 if + c 3 if 

~7 ~ (i^vT 2 ; 



(30) 



where Z is the compressibility factor, rj = p£, with £ = 
(7r/4)(xiCTj +X2CT2), is the total packing fraction, and the 
coefficients c\, C2, and C3 are obtained by identification 
with the corresponding coefficients which show up in the 
virial series. Specifically, in the present case one has 



B „ 
ci = - - 2, 



C n B 

C > = e- 2 T 



D C B 
(31) 

In Fig. [7] we present an illustrative comparison between 
the results for the compressibility factor of two binary 
NAHD mixtures as a function of the packing fraction 
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FIG. 5. Same as in Fig. [TJ but for A = 0.3. 

as derived from the RVE, Eq. (f30|) . and those obtained 
by MC simulation.— In both mixtures the size ratio is 
q = 0.7 and a negative nonadditivity A = —0.2 (with 
X\ = 0.5) and a positive nonadditivity A = 0.2 (with 
x\ = 0.4) have been considered. For comparison, the 
results stemming out of the compressibility factors cor- 
responding to the different theoretical approximations 
mentioned in section IIIII arc also included in this figure. 
As discussed in Ref. l35l for the actual calculations using 
the compressibility factors corresponding to the different 
theoretical approaches, one needs to specify the contact 
values of the one-component system for the Hamad and 
the SHY approaches and those of an additive hard-disk 
mixture in the MIX1 and mMIX theories. For the former 
we have used an accurate proposal by Luding ) 49 ' 50 while 
for the latter we have considered the quadratic approxi- 
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FIG. 6. Same as in Fig. [TJ but for A = 0.5. 



mation proposed in Ref. [5ll. complemented with Luding's 
one-component value.— 

It is clear that in the case of the mixture with negative 
nonadditivity, the best agreement is provided by both 
the RVE and the Hamad compressibility factor, followed 
by the SHY compressibility factor. In fact, the former 
two are hardly distinguishable. On the other hand, for 
positive nonadditivity it is the SHY compressibility fac- 
tor the one that provides the best agreement, followed 
by both the RVE and the MIX1 compressibility factor. 
These latter two are virtually indistinguishable. 
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FIG. 7. Plot of the compressibility factor Z versus the to- 
tal packing fraction r\ for NAHD mixtures with A = —0.2, 
ail<j\ = 0.7, x\ = 0.5 (top panel) and A = 0.2, (72/ci = 0.7, 
xi = 0.4 (bottom panel). The dotted (green) lines corre- 
spond to the original MIX1 theory, the dash-dot (pink) lines 
correspond to the mMIXl theory, the dashed (blue) lines cor- 
respond to Hamad's proposal, the thick solid (red) lines cor- 
respond to the SHY proposal, and the thin solid (black) lines 
correspond to the RVE, Eq. (f30|) . The symbols are our MC 
data.— Note that the RVE and the Hamad curves are prac- 
tically indistinguishable in the top panel, while the RVE and 
the MIX1 curves are practically indistinguishable in the bot- 
tom panel. 



VI. CONCLUDING REMARKS 

In this paper we have reported MC calculations of the 
fourth virial coefficients of asymmetric NAHD mixtures 
over a rather wide range of size ratios q and values of the 
nonadditivity parameter A. These results complement 
those reported earlier— for symmetric mixtures (q = 1) 
and, as illustrated in the case of the RVE and the mix- 
tures discussed in Sec. [Vj may prove useful for the devel- 
opment of new equations of state for NAHD mixtures. 
In particular, one could also consider using the avail- 
ability of the fourth virial coefficients provided in this 
paper to derive another approximation to the compress- 
ibility factor of asymmetric NAHD mixtures via the y- 
cxpansion proposed by Barboy and Gelbarti 52 ' 53 Here we 
have mainly used the data to assess the merits of differ- 



ent theoretical approaches leading to the thermodynamic 
properties of NAHD mixtures with respect to their per- 
formance in the prediction of the values of the fourth 
virial coefficients. 

One immediate conclusion is that none of the exist- 
ing theories can account for all the features observed in 
the MC data. In contrast with what happened in NAHS 
mixtures^ here the theoretical approach by Hamad^ 
outperforms all the rest. In this regard, it is somewhat 
striking that its very good performance concerning Z?m2 
and -D1222 is not also found for -D1122, where the SHY 
proposal does the best overall job. In any case, the com- 
parison we have presented is only indicative of the perfor- 
mance with respect to the fourth virial coefficients, but 
the full assessment will have to do with the compressibil- 
ity factor and with the issue of fluid-fluid demixing. We 
plan to address these points in the near future. 
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